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Abstract 


C  1 1-,.'*  V- 

We  investigate^the  linear  stability  of  plane  Couette  flow  of  an  upper  convected  Maxwell 
fluid  using  a  spectral  method  to  compute  the  eigenvalues.  No /instabilities  are  found. 
This  is  in  agreement  with  the  results  of  Ho  and  Denn^flpand  Leje  and  Finlayson  [2]'  but 
contradicts^proofs^of  instability  given  by  Gorodtsov  and  Leo^ov^[3]  and  Akbay  and 
Frischmann  [4,5]*  The  errors  in  those  arguments  are  pointed  out.^We  also  find  that  the 
numerical  discretization  can  generate  artificial  instabilities  (see  also  |l,6]*  The  qualitative 
behavior  of  the  eigenvalue  spectrum  as  well  as  the  artificial  instabilities  is  discussed. 
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SIGNIFICANCE  AND  EXPLANATION 


In  the  extrusion  of  polymers  from  a  pipe,  an  instability  known  as  melt  fracture  is 
frequently  observed.  There  is  considerable  controversy  concerning  the  question  whether 
this  instability  develops  in  the  inflow  region,  the  outflow  region  or  in  the  pipe  itself  and 
which  physical  effects  are  causing  it.  One  of  the  hypotheses  has  been  that  it  may  arise  from 
an  instability  in  parallel  shear  flow,  and  that  somehow  a  large  ratio  of  the  first  normal  stress 
difference  to  the  shear  stress  (or  a  similar  quantity)  is  responsible  for  causing  instability. 
In  this  paper,  the  authors  study  a  particular  model  of  a  viscoelastic  fluid,  and  the  results 
indicate  that  plane  Couette  flow  of  such  a  fluid  is  always  linearly  stable.  The  results 
therefore  question  the  explanation  of  melt  fracture  mentioned  above. 

The  paper  also  shows  that  artificial  instabilities  can  arise  from  numerical  discretiza¬ 
tion.  This  observation  is  relevant  in  the  context  of  numerical  simulations  of  flows  of 
viscoelastic  fluids,  which  have  encountered  many  problems. 
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LINEAR  STABILITY  OF  PLANE  COUETTE  FLOW 


OF  AN  UPPER  CONVECTED  MAXWELL  FLUID 

Michael  Renardv12  and  Yuriko  Renardy1 2 

1.  Introduction 

In  polymer  processing,  instabilities  known  as  melt  fracture  are  frequently  observed 
when  the  flow  rate  exceeds  a  critical  value  (see  [7,8]  for  reviews).  A  satisfactory  explanation 
for  these  instabilities  has  yet  to  be  found,  although  there  are  a  number  of  conjectures. 
There  may  be  more  than  one  mechanism  leading  to  these  instabilities.  For  example,  there 
is  some  discussion  whether  the  instability  originates  in  parallel  shear  flow  or  in  the  inflow 
or  outflow  region.  In  either  case,  it  is  not  understood  what  the  physical  mechanisms  are 
and  which  features  a  fluid  model  must  or  must  not  have  in  order  to  predict  the  instability. 

In  this  paper  we  study  the  linear  stability  of  plane  Couette  flow  of  an  upper  convected 
Maxwell  fluid.  In  our  study  we  assume  that  stability  can  be  determined  from  analyzing 
the  eigenvalues  of  the  differential  equation.  While  this  is  well  known  to  be  true  for  the 
Navier-Stokes  equations,  there  is  no  proof  of  this  for  the  equations  describing  the  upper 
convected  Maxwell  fluid.  For  a  discussion  of  possible  problems  that  can  arise  in  a  general 
context,  see  e.g.  section  4.4  of  [9]. 

Tlapa  and  Bernstein  [10]  have  shown  that  Squire’s  theorem  holds  for  the  upper  con¬ 
vected  Maxwell  fluid  (it  is  well  known  [11]  that  this  is  not  the  case  for  non-Newtonian  fluids 
in  general),  and  therefore  we  can  restrict  our  attention  to  two-dimensional  disturbances. 
There  is,  however,  a  special  class  of  disturbances  for  which  Squire’s  transformation  degen¬ 
erates,  and  which  therefore  requires  a  separate  investigation.  An  incorrect  investigation  of 
this  class  of  disturbances,  predicting  instability,  was  given  by  Akbay  and  Frischmann  [4,5]. 
We  shall  show  in  section  2  below  that  in  fact  these  disturbances  are  always  stable. 

Gorodtsov  and  Leonov  (3]  give  an  exact  analytical  solution  for  the  case  of  zero 
Reynolds  number.  They  find  that  in  addition  to  a  stable  continuous  spectrum,  there 
are  exactly  two  eigenvalues  for  each  value  of  a  (the  wave  number  in  the  streamwise  di¬ 
rection).  These  eigenvalues  are  also  stable.  Since  the  perturbation  introduced  by  a  small 
Reynolds  number  is  singular,  it  can  not  be  concluded  from  this  result  that  flows  at  small 
Reynolds  number  are  stable.  Gorodtsov  and  Leonov  give  an  incorrect  proof  of  instability 
at  finite  Reynolds  numbers  (see  the  discussion  in  section  2  below). 

Denn  and  his  coworkers  [1,6,12]  have  studied  the  stability  of  plane  Poiseuille  flow  of 
an  upper  convected  Maxwell  fluid  using  a  numerical  scheme  based  on  the  shooting  method. 
They  found  12]  that  viscoelastic  effects  are  destabilizing  at  high  Reynolds  number  and  the 

1  Sponsored  by  the  United  States  Army  under  Contract  No.  DAAG29-80-C-0041. 

2  Supported  by  the  National  Science  Foundation  under  Grant  Nos.  MCS-8215064  and 
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critical  Reynolds  number  is  decreased  by  a  small  amount  of  elasticity  (for  Couette  flow  of 
a  Newtonian  fluid,  there  is  no  critical  Reynolds  number  (13]).  At  low  Reynolds  numbers, 
no  instabilities  were  found  [l],  but  the  numerical  method  led  to  artificial  instabilities  [1,6]. 
The  recent  study  of  Lee  and  Finlayson  [2 ■ ,  which  uses  a  similar  numerical  method  to  study 
both  Poiseuille  and  Couette  flow,  confirms  the  absence  of  instabilities  at  low  Reynolds 
numbers. 

In  this  paper,  we  use  a  spectral  method  for  the  discretization  and  a  matrix  eigenvalue 
solver  for  computing  the  eigenvalues:  i.e.,  in  contrast  to  the  earlier  studies,  we  compute 
all  the  eigenvalues  of  the  problem,  rather  than  focussing  on  a  few  individual  eigenvalues. 
This  allows  some  insight  into  the  overall  qualitative  behavior  of  the  spectrum  as  well  as  the 
qualitative  nature  of  artificial  instabilities  produced  by  the  discretization.  Again,  no  insta¬ 
bilities  are  found.  In  fact,  the  eigenvalues  at  low  Reynolds  number  and  high  Weissenberg 
number  appear  to  be  rather  uninteresting,  although  they  are  difficult  to  compute.  The 
growth  rates  turn  out  to  be  close  to  the  constant  -1/2U''  (W  is  the  Weissenberg  number). 
At  high  Reynolds  number,  we  found  that  a  small  amount  of  elasticity  has  a  destabilizing 
effect,  but  we  did  not  find  instability. 

As  far  as  the  explanation  of  melt  fracture  is  concerned,  these  results  leave  three 
possible  conclusions: 

(a)  Melt  fracture  does  not  originate  in  parallel  shear  flow. 

(b)  Melt  fracture  is  a  finite  amplitude  effect  which  can  not  be  explained  by  linear  stability 
analysis. 

(c)  Melt  fracture  can  not  be  explained  by  the  upper  convected  Maxwell  model. 

In  connection  with  the  last  alternative,  it  may  be  relevant  that  many  models,  but  not 
the  upper  convected  Maxwell  model,  exhibit  instabilities  associated  with  a  change  of  type 
[14].  Also,  a  corrected  analysis  for  the  special  class  of  perturbations  studied  by  Akbay  and 
Frischmann  seems  to  show  instability  (not  associated  with  a  change  of  type)  for  the  lower 
convected  Maxwell  model  (U.  Akbay,  private  communication). 


2.  The  governing  equations  and  some  analytical  results 

The  flow  considered  is  between  two  infinite  parallel  plates  located  at  y  —  ±1,  which  are 
moving  in  the  i-direction  with  velocities  +1  and  -1.  The  equations,  given  in  dimensionless 
form,  are  the  equation  of  motion 

R\  ~  +  (u  •  V)uj  =  -  Vp4-  div  T, 

1  ot 


div  u  =  0, 

and  the  constitutive  law  for  the  upper  convected  Maxwell  fluid 

T~W  ~  +  (u  •  V)T  -  (Vu)T  -  T(Vu)r]  -  Vu  -  (Vu)T. 
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(1) 


•  y-.-\ 


(2) 


Plane  Couette  flow  is  the  trivial  solution 

/ 2W  10 
Uo  =  (y,0,0),  T0  =  1  0  0 

\  0  0  0 

To  this  we  add  a  small  perturbation 

u  -  u0  -  v(y)et'ar-*v,Jf+<r‘,  T  =  T0  -  S(y)e,'a*+i',*+‘T/, 


(3) 


(4) 


and  the  equations  are  linearized  with  respect  to  the  components  of  v  =  (u.v,w)  and  S. 
Tlapa  and  Bernstein  i  10  show  that  Squire’s  transformation 

(o  )2  -  o2  —  /32,  ,3  =0.  o'  la'  =  a/a,  a'  u'  =  au  -+  @w,  v  -  r,  w'  =0, 


a'W’  =  aW,  a' R*  =  aR 


(5) 


can  be  used  to  transform  the  problem  to  the  two-dimensional  case.  However,  this  breaks 
down  for  the  case  where  au  +  0w  —  0,  and  hence  this  case  requires  a  separate  investigation. 
If  au  4  (3w  =  0,  it  follows  from  the  incompressibility  condition 


iau  4  v'(y)  +  i/3vu  =  0, 


(6) 


and  the  boundary  conditions  at  y  =  ±1  that  v  -  0,  and  equations  (10)  and  (12)  in  (lOj 
lead  to 

RS(o  4  iay)w  =  w"  4  2 iaWw'  -  (a2  4  /?2  4  2a2W2)u>,  (7) 


where 

S  =  1  +  W(o  +  iay), 

and  an  identical  equation  for  u.  In  (7),  we  set  w  =  e~taWy<j),  and  obtain 

RS(o  +  iay)<t>  =  <t>"  -  (a2  4-  (32  4  a2W2)<f>. 


(8) 


(9) 


For  o  -  (3  —  0,  this  agrees  with  equation  (6a)  in  [5].  We  now  multiply  (9)  by  S<£,  integrate 
from  -1  to  1  and  integrate  by  parts  using  the  boundary  conditions  </>(l)  =  <£(-l)  =  0. 
This  yields 


*/: 


(a  -  iay)!5|2jd»j2  dy 


J  S(a 2  4  Q2  4  a2W2)\<t>:2  dy. 


4  S'<t>'4>  dy 


(10) 


Taking  real  parts,  we  obtain 

RReoj  ':S'2:<t>l2  dy  —  -(1  -r  W’  Re  a)  J  \<t>r2  -  (a2  +  $2  +  a2W2)|d>|2 

3 


dy 


*  S 


-Jl 


-aW  Im 


J  <t>  <t>  dy. 


Using  the  inequality 


aW \<}>'4>\  <  Ult'f  -  a2W2\<t>\2), 
£• 


we  find  that  the  right  hand  side  of  (11)  is  negative  as  long  as  Re  a  >  -l/(2W).  Hence  the 
left  hand  side  must  also  be  negative,  and  we  conclude  that  Re  a  must  be  negative. 

We  now  turn  to  the  study  of  two-dimensional  disturbances.  For  this  case,  we  introduce 
a  stream  function  by  ti  =  <t>' ,  v  =  —i<x(f>.  The  equation  for  the  stream  function  reads  as 
follows  (see  [3,12]): 

<t> ^  +  bz{y)<t>"'  +  bz  (y)<t>n  +  h(y)d>'  +  f>o(y)<p 

=  SR(yta  +  o)(4>"  -  a2d>),  (13) 


with  boundary  conditions 


<*(-!)  =  <*'(-!)  =  <*(!)  =  «*'(!)  =  0- 


Here  5  is  given  by  (8)  and 


b3(y)  =  2iaW 


S  -  1 


b2{y)  =  -2a2  -  2 a2W 


6,(9)  - 

6oW  =  a.  +  2^3^i ±  + 


An  equivalent  form  is 


\S2-j-=  -  2iaWS^~  -  2 a2W2  -  a2S2]  +  2iaW -  a2  -  2a2W2]d> 

1  dy2  dy  J  vdy2  dy 

-S3R(tay  +  <r)(^2  -  a2)0  =  0.  (16) 

The  fact  that  the  differential  operator  can  be  factored  for  R  =  0  was  used  by  Gorodtsov 
and  Leonov  [3)  to  obtain  an  exact  solution  for  this  case.  The  line  segment  from  -  pp  -  ia 
to  -  ^  +  ia  is  a  continuous  spectrum  arising  from  the  singular  character  of  the  equation 
when  S  =  0.  Apart  from  that,  there  are  exactly  two  eigenvalues  for  each  value  of  a,  for 
which  Gorodtsov  and  Leonov  give  an  exact  expression. 

Gorodtsov  and  Leonov  also  claim  to  find  eigenvalues  with  positive  real  parts  when  R  ■£ 
0.  Their  analysis  is  based  on  their  equation  (4.2)  which  they  derive  as  an  approximation: 

<t>(A)  -4-  2i aW<t>'"  -  2 a2(l  ±  (W2yy/ 2  -  (W2t)<t>"  -  2a2lT(ta  ±  (Wy/2)<t>' 


■ 


+  o4(  1  ±  2V2 tW2y  -  2 (W2^)<j>  =  0. 


(17) 

Their  instability  proof  is  based  on  finding  eigenvalues  £  with  positive  imaginary  parts  for 
this  equation.  In  fact,  however,  all  eigenvalues  for  (17)  are  real.  To  see  this,  we  multiply 
(17)  by  4>  and  integrate  from  -1  to  1.  After  some  integrations  by  parts,  we  obtain 


J  o"  2  -  2iaWQ"o'  -s-  2a2'<<t>'\2  -  2to3M’<s,<j  -  o4(l  ±  2\^2eW2y)\<p\2  dy 

?:2a2\Y2eV2  J  y\<S>"\2  dy  -  2a2W,2££  J  ti>' 2  +  a2  &2  dy  =  0.  (18) 

Since  all  other  terms  in  this  equation  are  real,  £  must  also  be  real. 

For  q  =  0.  equation  (13)  reduces  to 

<*>(4)  =  (1  +  Wo)o(t>".  (19) 


The  case  a  =  0  is  a  special  case  as  far  as  boundary  conditions  are  concerned,  because 
v  —  -ia<t>  -  0  does  not  imply  4>  =  0.  However,  the  boundary  condition  <f>  =  0  is  obtained 
if  the  perturbation  to  the  flow  rate  f]_  j  u(y)  dy  is  constrained  to  be  zero  and  a  uniform 
pressure  gradient  in  the  z-direction  is  allowed  (alternatively,  one  can  allow  a  nonzero  flow 
rate,  but  no  pressure  gradient;  in  that  case,  one  obtains  a  different  boundary  condition). 
Equation  (19)  with  boundary  conditions  (14)  has  the  eigenvalues 


a 


-1  ±  y/l+AWo'/R 
2  W 


(20) 


where  o'  =  (1  +  oW)o  is  given  by  -n2 zr2,  n  €  N  or  by  -/32  where  is  a  nonzero  root 
of  0  =  tan  0.  Except  for  a  finite  number,  all  roots  given  by  (20)  have  real  part  -  . 

Moreover,  the  convergence  of  the  real  parts  to  -  ^  is  most  rapid  when  W/R  is  large. 


3.  Numerical  Results 

For  the  numerical  solution,  we  multiplied  equation  (13)  by  S2,  and  discretized  the 
equation  and  boundary  conditions  by  the  Chebyshev  r-method  (see  (I5j).  Since  o  occurs 
in  powers  up  to  the  fourth  in  the  equation,  this  leads  to  a  matrix  eigenvalue  problem  of 
the  form 

det  (o*A4  +  +  o2 A2  +  o A\  +  Aq)  —  0.  (21) 

Since  four  equations  contain  the  boundary  conditions,  the  number  of  eigenvalues  is  4( N  -  4) 
when  A'  Chebyshev  modes  are  used.  We  rewrite  (21)  as  a  first  order  equation  with  a  matrix 
4  times  the  size  of  the  A,,  which  is  solved  using  the  NAG  routine  F02GJF.  This  procedure 
allows  us  to  compute  all  the  eigenvalues  rather  than  just  individual  ones.  It  is,  however, 
quite  sensitive  to  round-off  error,  and  calculations  at  high  values  of  aW  had  to  be  done  in 
quadruple  precision  (on  a  VAX  11/780).  In  a  few  cases,  we  obtained  improved  accuracy 


for  individual  eigenvalues  by  using  a  large  number  of  Chebyshev  modes  and  a  Newton 
iteration  based  directly  on  (21).  For  this,  quadruple  precision  is  not  necessary. 

The  program  was  tested  against  the  exact  result  given  by  (20).  We  also  compared  our 
results  with  those  of  Lee  and  Finlayson  !2i.  They  give  a  table  listing  eigenvalues  computed 
for  R  ~  0.25.  IV  -  1,  a  =  15  (in  comparing  their  results  with  ours  it  must  be  noted 
that  they  consider  Couette  flow  on  the  interval  0.1  rather  than  I  — 1,1,  and  that  their 
eigenvalues  are  wave  speeds  rather  than  growth  rates:  we  have  transformed  their  results 
to  our  present  notation).  The  eigenvalues  of  Lee  and  Finlayson  were  computed  using  an 
approximate  equation.  As  Lee  and  Finlayson  observe,  this  approximate  equation  produces 
the  imaginary  parts  of  the  eigenvalues  almost  perfectly,  but  approximates  the  real  parts 
rather  poorly,  leading  to  artificial  instabilities  at  higher  Weissenberg  numbers.  This  is 
also  reflected  in  the  following  comparison  of  our  results  (using  80  Chebyshev  modes)  with 
theirs. 

Eigenvalues  with  lowest  imaginary  parts  Eigenvalues  given  in  table  VI  of  [2] 
computed  by  our  program 


0.85349  -  14.646? 
0.50964  -  33.021? 
0.50824  *  36.766i 
0.50694  -  39.783? 
0.50579  ±  42.414? 
0.50480  x  44.793? 
0.50398  ±  46.986? 
0.50332  *  49.038? 
0.50280  ±  50.975? 
0.50253  i  52.824? 
0.50283  .t  54.631? 
0.50387  ~  56.480? 
0.50518  i  58.442? 
0.50623  2  60.537? 
0.50683  2  62.755? 
0.50700  r  65.074? 
0.50685  z  67.479? 
0.50646  r  69.956? 
0.50595  r  72.494? 
0.50536  r  75.085? 
0.50478  »  77.720? 
0.50421  r  80.395? 
0.50369  *  83.105? 
0.50322  2.  85.846? 
0.50281  2  88.615? 
0.50245  2  91.408? 
0.50214  2  94.224? 
0.50187  2  97.060? 


-0.960  ±  33.021? 


-0.834  +  44.796? 


-0.543  -  62.757? 


-0.432  -  75.084? 


-0.417  -  91.407? 
-0.420+  94.224? 


0.50164 

0.50144 

-0.50128 

-0.50119 

-0.50135 

0.50236 

-0.50579 

0.51427 


99.915? 

102.79? 

105.67? 

108.58? 

111.49? 

114.42? 

117.37/ 

120.35? 


-0.429+  105.67* 


-0.441  a:  120.30? 


The  sparseness  of  the  right  column  shows  that  the  list  of  eigenvalues  given  by  Lee  and 
Finlayson  is  far  from  complete.  As  can  be  seen,  the  real  parts  of  all  our  eigenvalues  fall 
almost  exactly  on  -  ,  except  for  the  first  pair.  This  one  pair  of  eigenvalues  behaves  in  a 

qualitatively  different  fashion  from  the  others.  It  corresponds  to  the  two  eigenvalues  found 
by  Gorodtsov  and  Leonov  :3‘  for  R  -  0.  The  occurrence  of  the  eigenvalues  in  complex 
conjugate  pairs  is  due  to  the  reflection  symmetry  of  the  Couette  flow  problem  about  the 
origin.  For  imaginary  parts  beyond  120.  the  numerical  accuracy  of  our  results  deteriorates, 
as  already  evident  in  the  deviation  of  the  real  parts  of  the  last  few  eigenvalues  in  the  table 
from  -  jpy  •  We  recomputed  the  last  eigenvalues  in  the  table  above  with  120  Chebyshev 
polynomials  and  found  -0.50074  ±  120.30?. 

We  did  a  number  of  calculations  at  R  —  1,  a  =  1  and  various  Weissenberg  numbers. 
The  following  table  shows  the  eigenvalues  with  the  smallest  imaginary  parts. 


W  =  0.2 


W  =  2 


W  =  20 


W  =  50 


4.9687  ± 
2.5311  ± 
2.4949  ± 
2.5063  ± 
2.4971  r 
2.5026  ± 
2.4984  * 
2.5014  ± 
2.4990  ± 
2.5009  2 


0.6462? 

6.3370? 

9.8315? 

13.726? 

17.155? 

20.861? 

24.300? 

27.940? 

31.389? 

34.997? 


-0.36985 

0.34582 

-0.27904 

-0.24678 

-0.25892 

-0.24631 

-0.25396 

0.24762 

0.25226 

0.24841 


±  0.7740? 
±  2.1862?' 
±3.5116? 
±  4.5354? 
±  5.6822? 
±  6.7446? 
±  7.8720? 
±  8.9500? 
±  10.073? 
±  11.160? 


-0.026286 

-0.025106 

-0.025151 

-0.025212 

-0.025294 

-0.025408 

-0.025527 

-0.025762 

-0.026085 

-0.026198 


±  0.9750?' 
±  3.8928? 
±  4.1968? 
±  4.4433? 
±  4.6597? 
±  4.8565? 
±  5.0397? 
±  5.2134? 
±  5.3910? 
±  5.5778? 


-0.010207  ± 
-0.010003  ± 
-0.010004  ± 
-0.010004  ± 
-0.010005  ± 
-0.010005 ± 
-0.010006 ± 
-0.010006  ± 
0.010007  ± 
-0.010008  ± 


0.9900? 

6.3345? 

6.5290? 

6.6874? 

6.8269? 

6.9540? 

7.0720? 

7.1830? 

7.2884? 

7.3891? 


The  first  pair  of  eigenvalues  corresponds  to  those  of  Gorodtsov  and  Leonov.  They 
approach  the  continuous  spectrum  (i.e.  the  line  segment  from  - ^  -  ?a  to  +  to)  as 
W  —  0.  For  large  W .  the  real  parts  approach  -  ^  and  the  imaginary  parts  tend  to  ±a. 
For  R  =  0.  this  behavior  can  be  deduced  from  the  formula  given  in  (3]. 

The  remaining  eigenvalues  are  basically  lined  up  along  the  line  Re  a  —  —  ^ .  We  see 
in  the  first  two  columns  that  the  real  parts  approach  —  ^  as  the  imaginary  parts  increase, 
this  is  in  fact  suggested  by  the  analysis  for  a  =  0.  The  eigenvalues  in  the  third  column 
appear  to  move  away  from  -  ^ ,  but  this  trend  is  reversed  as  the  imaginary  parts  increase 
further.  For  W  ~  50,  a  large  number  of  Chebyshev  modes  (we  used  120)  was  required 
to  obtain  sufficiently  good  results.  All  the  eigenvalues  at  W  =  20  and  in  particular  at 
W  -  50  have  real  parts  close  to  -  An.  This  suggests  that  the  real  parts  of  the  eigenvalues 


approach  -  ^  at  large  W.  as  is  indeed  the  case  for  a  =  0.  The  results  tabulated  above 
for  a  -  15  suggest  that  the  real  parts  of  the  eigenvalues  also  tend  to  -  ^  at  large  a.  The 
spacing  between  the  imaginary  parts  of  the  eigenvalues  decreases  with  increasing  W ,  while 
the  imaginary  part  of  the  second  eigenvalue  first  decreases  and  then  increases  again.  This 
behaviour  would  be  expected  from  an  analysis  of  characteristic  wave  speeds  [14]. 

While  the  actual  spectrum  of  the  differential  equation  shows  no  instabilities,  the  nu¬ 
merical  approximation  does.  In  the  numerical  calculations  at  low  Reynolds  numbers,  we 
can  distinguish  four  clearly  separated  sets  of  eigenvalues: 

1.  Spurious  eigenvalues:  By  this  we  mean  eigenvalues  of  the  discretized  problem  which 
do  not  approximate  those  of  the  differential  equation  even  qualitatively  and  lie  in  a 
totally  different  part  of  the  complex  plane. 

2.  The  two  “Gorodtsov-Leonov”  eigenvalues. 

3.  Eigenvalues  approximating  the  remainder  of  the  discrete  spectrum. 

4.  Eigenvalues  approximating  the  continuous  spectrum,  i.e.  the  line  segment  from  —  ^  - 
ia  to  —  vv  “b  Ia- 

In  all  our  calculations  we  found  four  spurious  modes.  These  spurious  modes  exist  even 
in  the  Newtonian  case  [15]  and  they  can  be  unstable.  A  certain  number  of  the  eigenvalues 
in  the  third  group  gives  good  approximations  to  those  of  the  differential  equation,  but  as 
the  imaginary  parts  increase,  the  accuracy  ultimately  deteriorates.  Since  the  real  parts 
are  small  compared  to  the  imaginary  parts,  they  are  particularly  affected  and  they  can 
in  fact  deteriorate  to  the  point  where  they  have  the  wrong  sign,  thus  creating  artificial 
instabilities.  Not  surprisingly,  this  is  most  likely  to  happen  at  high  Weissenberg  numbers. 
The  approximation  to  the  continuous  spectrum  is  generally  poor.  This  is  not  surprising. 
The  method  used  here  approximates  isolated  eigenvalues  with  C°°-eigenfunctions  with 
infinite  order  accuracy  [15],  but  this  is  not  the  case  for  a  continuous  spectrum.  The 
approximation  to  the  continuous  spectrum  is  best  near  the  ends  and  worst  near  the  middle. 
Again  artificial  instabilities  develop  at  high  Weissenberg  numbers.  We  believe  that  the 
“spurious”  modes  of  Ho  and  Denn  [l]  also  result  from  poor  approximation  to  the  continuous 
spectrum.  The  following  table  shows  the  number  of  numerically  unstable  eigenvalues  as  a 
function  of  the  number  of  Chebyshev  modes  at  R  =  1,  a  =  1  and  W  —  20. 


No.  of  modes  30 

40 

50 

60 

80 

Group  1  2 

2 

2 

2 

2 

Group  3  8 

10 

14 

16 

22 

Group  4  17 

19 

19 

13 

0 

We  see  that  the  instabilities  from  the  continuous  spectrum  disappear  if  a  sufficient 
number  of  Chebyshev  modes  is  used,  while  those  from  the  discrete  spectrum  only  get 
shifted  towards  higher  imaginary  parts,  and  the  number  of  unstable  eigenvalues  actually 
increases. 

The  equations  for  steady  flow  of  an  upper  convected  Maxwell  fluid  undergo  a  change 
of  type  when  RW  —  1  -t  IV2  j  1 4 ] .  All  the  results  reported  above  are  subcritical,  i.e. 
RW  <  1  +  W2.  As  an  example  of  a  subcritical  case,  we  looked  at  R  =  10,  W  =  2,  a  =  5. 
The  eigenvalues  with  the  lowest  imaginary  parts  are  as  follows: 


-  0.27302  ±  0.1768? 

-0.28615  X  0.8463? 

-0.25010  ±  0.9283? 

-0.24455  ±  1.5281? 

0.25267  x  2.0821? 

0.25051  -  2.5991  / 

-  0.24896  x  3.0702? 

-0.25073  x  3.5151? 

-0.24993  x  3.9371/ 

0.24991  4.3361/ 

-0.25015  x  4.7198/ 

-0.35520  x  4.7784/ 

0.24969  x  5.0858? 

0.25060  ..  5.4386? 

The  change  of  type  means  that  in  a  part  of  the  domain  the  speed  of  the  fluid  exceeds 
a  characteristic  wave  speed.  This  manifests  itself  in  the  imaginary  parts  of  the  eigenvalues, 
which  now  “fill  out”  the  whole  real  axis  including  the  interval  (-a,  a),  rather  than  being 
separated  by  a  gap  in  the  middle.  The  qualitative  behavior  of  the  real  parts,  however,  is 
unchanged;  they  are  still  close  to  except  for  the  twelveth  pair,  which  represents  the 

Gorodtsov-Leonov  eigenvalues. 

We  did  a  few  calculations  at  a  high  Reynolds  number  ( R  =  10000,  a  =  1).  In  the 
Newtonian  case,  the  following  asymptotic  formula  for  the  least  stable  eigenvalues  holds  for 
R  >  oo  (see  [16],  section  31.1): 

<7  =  -1.0626a2/37r1/3it(a-  4.1288a2/3R~1/3).  (22) 

At  R  =  10000,  o=l,  this  is  equal  to  -0.04932  X  0.80836?.  The  following  table  shows  our 
computed  results  at  various  values  of  W . 

W  =  0.01  W  =  0.1  W  =  0.5  W  =  2 

-0.05205  x  0.81213?  -0.05167X0.81164?  -0.04984  X  0.80960?  -0.04157  X0.80409? 

We  see  that  elasticity  has  a  destabilizing  effect,  but  it  does  not  seem  to  lead  to  insta¬ 
bility.  Apart  from  the  spurious  modes,  we  found  no  artificial  instabilities  of  the  numerical 
discretization,  unless  very  few  Chebyshev  modes  are  used.  Even  though  the  numerical 
method  is  slower  to  converge  at  higher  Reynolds  numbers,  it  is  less  likely  to  produce 
unstable  eigenvalues. 
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